Qubit vitrification and entanglement criticality on a quantum simulator

Many elusive quantum phenomena emerge from a quantum system interacting with its classical environment. Quantum simulators enable us to program this interaction by using measurement operations. Measurements generally remove part of the entanglement built between the qubits in a simulator. While in simple cases entanglement may disappear at a constant rate as we measure qubits one by one, the evolution of entanglement under measurements for a given class of quantum states is generally unknown. We show that consecutive measurements of qubits in a simulator can lead to criticality, separating two phases of entanglement. Using up to 48 qubits, we prepare an entangled superposition of ground states to a classical spin model. Progressively measuring the qubits drives the simulator through an observable vitrification point and into a spin glass phase of entanglement. Our findings suggest coupling to a classical environment may drive critical phenomena in more general quantum states.

The Born rule, which states that the outcome of a measurement performed on a quantum state is a random variable whose probability distribution is determined by quantum theory, governs information transfer from a quantum system to its classical environment. While the Born rule is at play constantly all around us since all matter is fundamentally quantum, its effects are only evidenced in bulk due to the astronomical number of measurement events that occur at macroscopic lengths and time scales. In contrast, a quantum simulator, a programmable array of qubits, is an otherwise isolated quantum system that we can couple to its environment at will with tailor-made measurements of all or some of the qubits. As such, it allows us to study in detail how the quantum characteristics of a system change as we progressively measure its components.
If we were to pick the state of an ideal quantum simulator uniformly at random from the set of all states accessible to it, we would get a volume-law state: a state whose entropy of entanglement of an extensive subsystem, measured in bits, is proportional to the number of qubits. Measuring a single qubit in such a state removes roughly one bit of entanglement. One therefore expects that the entanglement of a volume-law state should decrease as a simple linear function of the number of measured qubits, yielding a classical unentangled state once we measure all the qubits. However, this is not always the case: the behaviour of entanglement in a quantum simulator can change dramatically and abruptly as we progressively measure its qubits, exhibiting the phenomenon of criticality 1 . Critical behaviour indicates a transition between two distinct phases of entanglement. Previous work revealed entanglement phase transitions in different settings, namely, in random ensembles of quantum circuits in which qubits undergo measurement at a finite rate [2][3][4][5][6][7][8] and in models with topological order [9][10][11] .
Detecting and characterizing such critical behaviour in experiments is challenging though. First, current quantum processors are faulty, limiting the quantum gates we can reliably implement and the number of error-free measurements we can obtain. Second, measuring the entanglement entropy of an arbitrary quantum state is hard. The straightforward approach requires tomographic reconstruction of the quantum state, which is intractable for quantum systems with many components. Finally, although theoretical models for entanglement phase transitions have been introduced in the context of random circuit ensembles 2,3,5,8 , these are only approximate in the experimentally relevant limit.
Here, we program two entanglement phases and the criticality between them on a quantum simulator of up to 48 superconducting qubits. We do this by implementing an ensemble of quantum circuits 1 that allow us to reliably generate volume-law states whose entanglement we can deliberately decrease with qubit measurements, experimentally determine the entanglement entropy, and capture the exact dependence of entanglement on measurements using a physical theory 1 (see "Measuring entanglement"). The pertinent physical theory is that of spin vitrification, i.e., the transition to a spin glass phase 12,13 . We detect the vitrification point, which agrees with spin glass theory. Our work shows measurements alone can trigger entanglement criticality, suggesting a classical environment could induce critical behaviour in more general quantum states.

Theory and model
Our experimental system is an array of superconducting qubits on a quantum simulator. To drive these qubits through an entanglement phase transition, we first execute quantum circuits that prepare highly entangled states. The circuits and theory that follow come from previous theoretical work 1 on entanglement phase transitions.
Each of our circuits implements a system of R linear equations on L Boolean variables using R + L = N qubits. (In practice, we use fewer qubits to implement our system on hardware. See "Circuit optimization" in the "Methods".) We can write the system as the matrix equation where B is an R × L Boolean matrix whose rows represent equations and columns represent variables, as shown in Fig. 1a. If B ij = 1, then equation i involves variable j. Otherwise, the entry is zero. By setting the elements of B according to some distribution, we get an ensemble of matrices. Each element of the parity vector y 2 0, 1 f g R fixes the parity of an equation and each x 2 0, 1 f g L is a solution to the system for a given y.
To implement the system of Eq. (1) using a quantum circuit, we organize the qubits in the simulator in two registers, as sketched in Fig. 1b,c: a "variable" register consisting of L variable qubits and a "parity" register consisting of R parity qubits. We input variable qubits in the |+i state and parity qubits in the |0i state into a circuit from the ensemble defined above. The initial state is therefore |ψi in = |+ i L |0i R . The state |ψi out at the output of the circuit is an entangled equal superposition of solutions x for each possible y (see "Quantum state and entanglement entropy" in the "Methods" for the exact state). For each y, the set of solutions x f g is unique. Each parity qubit at the output holds the parity of the variables that appear in the corresponding row of B (see Fig. 1b). Because these variable qubits are in a superposition of classical states, so is the parity qubit: it is 0 or 1 depending on the state of the variable qubits. The variable qubits are thus entangled with the parity qubit and contribute one bit of entanglement.
What is the total entanglement between variable and parity qubits? We quantify this with the entanglement entropy S, which counts the number of bits of entanglement between two parts of a quantum system. To answer our question above, we need to know how many possible vectors y there are in the superposition. For each of the rank(B) linearly independent rows of B, the corresponding component in y can be set to zero or one freely. There are then N y = 2 rank(B) possible vectors y which give solutions to Eq. (1). This means the entanglement entropy between the variable and parity qubits is Sð|ψi out Þ ∼ log 2 N y = rankðBÞ.

Entanglement phase transition
Each time we measure a qubit in a generic volume-law state, the system loses one bit of entanglement. This happens with every measurement, reducing the number of superposed configurations until just a single classical state remains.
The states |ψi out behave in a manifestly different manner. To prove this, we start by compiling enough linear equations into our circuits such that an output state is a superposition of N y = 2 rank(B) = 2 L vectors y. Then, we calculate the entanglement entropy after measuring a subset M of the parity qubits in our system in the computational basis. We label the measurement outcome y out,M and the partially measured state |ψ out,M . The state |ψ out,M still contains an equal superposition of solutions to Eq. (1), but the elements of the parity vector y that correspond to the subset M are fixed to the measurement outcome y out,M . For the same reasoning as in the previous section, there are N y out,M = 2 rankðB M Þ equally probable measurement The first linear equation in B compiled as a circuit. The gate labelled "1" is a CNOT followed by a SWAP, as shown in the inset of (c). We use this circuit design to construct the full quantum circuit for Bx = y. c The quantum circuit built from B. B ij = 1 corresponds to the two-qubit gate shown in the inset, whereas B ij = 0 corresponds to a SWAP gate. The output state |ψi out holds all y that yield solutions to Eq.
(1) and the corresponding x for the matrix B. Throughout all panels, the variable labels are in light blue and the parity/equation labels are in dark blue. outcomes y out,M , where B M are the rows of B that correspond to the M parity qubits. (In Fig. 1c, measuring the first three parity qubits would mean B M is the first three rows of B.) Since measuring the output state determines y out,M , there are N y =N y out,M = 2 LÀrankðB M Þ vectors y remaining in the state |ψi out,M . The entanglement entropy is then Sð|ψi out,M Þ ∼ log 2 ðN y =N y out,M Þ = L À rankðB M Þ. Therefore, the evolution of the entanglement entropy is given by rank(B M ) as a function of M. To obtain a size-independent control parameter, we define α ≡ |M|/L, the ratio of measured parity qubits to variable qubits.
A unique feature of states |ψ out is that we can obtain an exact result for the evolution of the entanglement entropy under measurements through a mapping to a classical spin model 1 . In this description, |ψ out is an entangled superposition of ground states x to a spin Hamiltonian with couplings defined by y (see "Methods"). We can then use the characteristics of the spin model both to predict the behaviour of the entanglement and, more importantly, to measure it on a quantum simulator, as we discuss in the next section.
To get concrete predictions for our experiments, we now specify the distribution we sample to populate the matrix B and define the ensemble of states |ψ out . We pick three distinct variables uniformly at random for each equation in Eq. (1), and we ensure there are no repeated equations. With this choice, we get an exact correspondence between our output state and the ground states to an instance of the unfrustrated 3-spin model, both of which are given by the solutions to Eq. (1) 14,15 . This model exhibits a phase transition at α c ≈ 0.918. For α < α c , our system corresponds to a paramagnet in the 3-spin model. We thus get a paramagnetic phase of entanglement. Here, rank(B M ) = |M| = Lα. This happens because there are few rows in B M , which makes it highly probable that they are linearly independent. The entanglement entropy of the quantum system after measuring |M| parity qubits scales linearly in both L and α: Sð|ψ out,M Þ ∼ L 1 À α ð Þ, i.e., the output state obeys a volume law. For α > α c , the system enters a spin glass phase. The qubits vitrify, turning into an entangled superposition of spin glass ground states. We thus get a glassy phase of entanglement. Now, rank(B M ) < |M| because there are many rows in B M and linear independence is lost. The entanglement entropy still scales linearly in L but decreases slower than linearly with increasing α.
We sketch how the measurement of parity qubits collapses the output state in Fig. 2. The two entanglement phases give rise to different behaviours. In the paramagnetic phase, measuring a parity qubit collapses its state to one of two equally probable values (Fig. 2b). This occurs when a parity qubit is independent of the other measured parity qubits, which is the case when B M has full rank. Measuring the parity qubit halves the number of possible vectors y remaining in the superposition, so the system loses one bit of entanglement entropy. In the spin glass phase, there is a finite probability that a parity qubit has a definite value before measurement (Fig. 2c). This occurs when previous measurement outcomes determine the measurement outcome for the next parity qubit, which begins when rank(B M ) < |M|. In this case, measuring the parity qubit does not change the number of possible vectors y, so the entanglement entropy remains the same.

Measuring entanglement
The exact correspondence between spin glass physics and entanglement established above and in ref. [1] gives us an efficient way to detect the entanglement phases and criticality on a quantum simulator. In our setup, the entanglement entropy is characterized by the spin glass order parameter [16][17][18] where 〈…〉 is an average over all the solutions x for Eq. (1) with matrix B M and a given parity vector y out,M , and x i is the i-th variable in x. We derive the identity that links this order parameter to the entanglement entropy in the "Methods". Therefore, while in most quantum systems quantifying entanglement is intractable, here we have direct access to the entanglement entropy through the spin glass order parameter, which is efficiently measurable (see "Methods").
The order parameter describes the tendency of the solutions x to take the same value on each variable. It is zero in the paramagnetic phase, which means a given variable does not have correlations across solutions. The outcome of the measurement of a parity qubit is independent of previous parity measurements in this phase. At α = α c , the variables abruptly become correlated across solutions, and the order parameter jumps to a finite value, eventually saturating to one. This implies solutions of the system are almost identical, differing on only a few variables. The outcome of the measurement of a parity qubit now depends on previous parity measurements.
In the language of physics, each variable can be thought of as one of L spins in a many-body system with |M| interactions. Then, each basis state in |ψ out,M of the variable qubits represents a spin configuration. In the paramagnetic phase where there are few interactions, these configurations have no correlation, leading to no order (q = 0). However, in the spin glass phase where |M| > Lα c , the interactions induce correlations across the configurations, leading to spin glass order (q > 0).
Using arrays of up to 48 qubits, our experimental results (Fig. 3) clearly reveal the two entanglement phases and the transition between them, and are in agreement with theory. The transition at a critical value α c becomes more abrupt with increasing system size, exactly as spin glass physics dictates. Finite-size scaling (see "Methods") of  (c). In the paramagnetic phase, measurement collapses the parity qubit, decreasing the entanglement entropy by one bit. In the spin glass phase, a finite fraction of the parity qubits are already collapsed from previous measurements, so measuring them does not affect the entanglement. An abrupt change between the two behaviours occurs at |M| = Lα c . the data using the scaling form qðαÞ = f ððα À α c, exp ÞL 1=ν exp Þ yields the experimental values for the critical measurement ratio α c, exp = 0:95 ± 0:06, which agrees with the theoretical value α c ≈ 0.918 15 , and the critical exponent ν exp = 2:5 ± 0:5.
We note that while there are some similarities between the spin glass order we find and those in other work on entanglement criticality [19][20][21] , there are a few differences. First, previous work focuses on states and circuits that respect certain symmetries, which play a role in creating a spin glass phase. In contrast, we do not need to impose symmetries. Second, there is an exact relation between the spin glass order parameter in our work and the entanglement entropy, which is not there in other works. Third, previous work focuses on spin glass order as a steady state property of the system, whereas our system goes into a spin glass state immediately after applying our circuit and measuring. Finally, as our results demonstrate, we can observe this spin glass order on existing quantum hardware.

Discussion
Since entanglement is a key resource for quantum computation, precise predictions and experimental verification of its possible behaviours in quantum devices under measurement are soughtafter. Our findings demonstrate that partial measurements of quantum states can alone give rise to intricate phenomena related to entanglement. Measurements can force qubits to vitrify, and hence realize the celebrated 13 spin glass phase of matter inside a quantum processor.
The spin glass quantum states implemented here are a subset of stabilizer states, an important class of states for quantum computation. Moreover, we already know that spin glass entanglement criticality is also present in more general classes of states than the ones studied here 1 . It is interesting to ask whether similar physics applies to entanglement in monitored quantum systems at large, giving rise to different types of nonanalyticity for the entanglement entropy.

Spin Hamiltonian
The output of our circuits provides x and y from Eq. (1), which we can relate to the ground states and couplings of a p-spin model 22 Fig. 1a). The indices of the nonzero elements in each row a ∈ B correspond to the spins which are part of an interaction. The Hamiltonian is: where a 1 , a 2 , a 3 refer to three distinct spins (the indices of the nonzero elements in a), J a are the couplings of the interaction vector J 2 ± 1 f g R , and the spins σ i form the spin vector σ 2 ± 1 f g L . The ground-state energy for this Hamiltonian is zero, which corresponds to J a σ a 1 σ a 2 σ a 3 = 1 for all a.
Using the mapping J a = ðÀ1Þ y a and σ i = ðÀ1Þ x i , we see that Eq. (3) is zero whenever y a = x a 1 + x a 2 + x a 3 mod 2 for all a, which is Eq. (1). This lets us express the number of ground states N GS to Eq. (3) in terms of y and B. Because each of the N y = 2 rank(B) vectors y has N GS ground states (out of a possible 2 L ), we find N GS = 2 LÀrankðBÞ . The ground-state entropy is S GS ðBÞ log N GS = L À rank ðBÞ Â Ã log2 (we take the natural logarithm).

Quantum state and entanglement entropy
After applying the circuit given by B (Fig. 1c) to our input state |ψ in , we get the following state (see ref. 1 for more details): This is a superposition of solutions x f g for each of the N y possible y. We then measure the state of the first |M| parity qubits to be y out,M . The resulting state is The entanglement entropy coincides with S GS (B M )-the groundstate entropy of the classical spin model-when we choose our initial B to satisfy rank(B) = L. In the limit of large system sizes, the expression for the averaged entropy density 23 is: where q(α) is the spin glass order parameter after performing the ensemble average using Eq. (2). Equation (7) establishes the exact correspondence between the entanglement entropy and the spin glass order parameter.

Quantum hardware
We used the ibm_washington, ibmq_brooklyn and ibm_hanoi quantum processors 24 for our experiments. We set the repetition delay to 0.00025s. We chose connected lines of qubits to take advantage of the one-dimensional structure of our circuits, while also having low measurement readout error and CNOT error rates at the time of job submission.

Circuit optimization
Because errors dominate the output in current quantum processors, we optimize our circuits to use as few gates as possible. As the CNOT is the native two-qubit gate on the IBM Q processors, we use the CNOT count N CNOT as our metric (we ignore the L single-qubit Hadamard gates we always need). We build our circuits using the matrix B M since it generates the solutions we use in Eq. (2) and requires fewer gates to implement than B. We then decompose SWAP gates in our circuits as SWAPði, jÞ = CNOTði, jÞ × CNOTðj, iÞ × CNOT ði, jÞ, ð8Þ where i and j are the qubits participating in the gate and CNOT(i, j) means qubit i controls the target qubit j. For the "1" gate in Fig. 1, using Eq. (8) reveals two consecutive CNOT gates with the same control and target, which we remove because they have no overall effect. As such, a one in the matrix requires two CNOTs while a zero requires three. How many qubits and CNOT gates do we need to build circuits such as in Fig. 1c using B M ? There are |M| = Lα linear equations, so we require N = L + |M| = L(1 + α) qubits. To calculate N CNOT , note that each row of B M contains p ones and L − p zeros. There are then 2p + 3(L − p) CNOTs per row of B M , where p = 3 for our model. Summing the CNOT count over all rows, we find N CNOT ðB M Þ = |M| 2p + 3ðL À pÞ ½ = 3L L À 1 ð Þα. By transforming B M using matrix row operations, we can reduce N and N CNOT . We begin by putting B M into row echelon form. Then for each row of the matrix (starting from the second), we find the index of the leading one, and add this row to the rows above it which have a zero at that index. These row additions generate more ones in the matrix, which we prefer because they require less gates than zeros to implement. We call the resulting matrix B M 0 . Note that solutions to Eq. (1) for a given parity vector remain unchanged under row operations.
where the omitted entries are zeros. The form of B M 0 helps us save qubits and gates. First, the matrix has size |M 0 | × L, with |M 0 | = rankðB M Þ. The resulting circuit requires N = L + rank(B M ) ≤ L + |M| qubits, fewer than the circuits built from B M when rank(B M ) < |M|. Second, notice that in Fig. 1c, there are gates for each entry of the matrix. By interspersing the parity and variable qubits instead of separating them, we can avoid including gates for the zeros to the left of the leading ones in B M 0 .
Each gate in the primitive circuit (Fig. 1b) exchanges the positions of the qubits it acts upon. The result is that a parity qubit exchanges positions with every variable qubit. However, only "1" gates contribute to the parity we want to measure. Therefore, once a parity qubit encounters all the "1" gates for its linear equation, we measure it in that position. We take advantage of this by altering the primitive circuit: we start the parity qubit at the top, reverse the gate sequence, and invert the control and target of the CNOTs in each "1" gate. Then, the locations of the leading ones in B M0 provide the end positions for measuring the parity qubits. For example, the parity qubit for the first row in Eq. (10) exchanges positions with all variable qubits before we measure it. The parity qubit for the second row exchanges positions with variable qubits 6, 5, 4, 3, and 2 before measuring, and so on.
We provide an upper bound for N CNOT ðB M0 Þ. We count the entries in B M0 to the right of (and including) the main diagonal. We assume the leading ones are all part of the main diagonal. With this assumption, the number of entries to the right of (and including) the leading ones in a P × Q row echelon matrix is: In our case, P = rank(B M ) and Q = L. There is at least a one per row (else the row would not be a part of B M0 ). Since zeros contribute more to N CNOT , we take the worst-case scenario where all other entries are zero. This implies there are P ones in the matrix, so there are Z = U(P, Q) − P zeros. Using these results, we get the upper bound: where we get the final inequality by maximizing the previous expression with respect to rank(B M ). Finally, rather than putting a variable qubit in its initial superposition | + i = H|0i as an input to the circuit, we apply the Hadamard gate H only when the corresponding variable first participates in a linear equation. (For example, in Eq. (10) variable 6 is first part of an equation in row 3.) Doing so reduces errors from trying to maintain superpositions for too long in current quantum processors. It also simplifies any SWAP gate involving a variable qubit in the state |0i. If we have qubits i and j with the latter in the state |0i, Eq. (8) reduces to SWAPði, jÞ| j in |0i = CNOTði, jÞ × CNOT ðj, iÞ: Our circuit optimization provides a significant savings compared to the circuits built from B M , which require N CNOT (B M ) = 3L(L − 1)α gates and N = L 1 + α ð Þqubits. In practice, our largest experiments (L = 24 and α ≳ 1) required an average of N CNOT ðB M 0 Þ≈600 gates, which is much less than the N CNOT (B M ) ≳ 1600 gates we would need if we used the B M matrices instead.

Error mitigation and shot count
In our experiments, we only keep measurement results (shots) x and y out,M 0 which satisfy B M 0 x = y out,M 0 . This provides significant error mitigation (see Supplementary Fig. 1) as we increase the number of qubits. For L = 8, 16,24 f g , we took 10000, 25000, 750000 f g shots per sample (see next section) to obtain our data. A. Fix a reference solution z that maps solutions from the parity vector y out,M 0 to the parity vector 0. We chose our reference to be the solution to B M 0 z = y out,M 0 whose binary form represents the smallest integer. Note that finding a reference is efficient. B. Map the saved solutions x associated with y out,M0 to x 0 = x + z. Now, B M 0 x 0 = 0. Remove any duplicates. Call this set X = x 0 f g. iii. Compute qðB M 0 Þ in Eq. (2) by uniformly sampling min 24,jX j ð Þ solutions from X, where we determined the number 24 yields a reasonable compromise between accuracy and quantum runtime. (b) Compute q(α) by averaging over qðB M 0 Þ for all B M 0 .
Following the procedure for each L produces the curves in Fig. 3. To calculate the order parameter in step iii), we want as many solutions x 0 as possible, but a finite number works, making the order parameter efficient to measure. For the classical simulation of q (the dashed lines in Fig. 3), we uniformly sample min 24, N GS À Á solutions to the equation B M x = 0, where N GS = 2 LÀrankðB M Þ is the total number of solutions. We did this by sampling random linear combinations of the basis vectors forming the null space of B M . This is also efficient.
We note that the small size of the sample X leads to appreciable artefacts in Fig. 3, such as a deviation of the order parameter from the expected value of zero at small α. Concomitantly, we notice the onset of finite-size effects at values of L and α for which min 24, N GS À Á ≈N GS . The dip of q at small α for L = 8 is such an effect. We nevertheless notice that, for all L and α, theory and experiment match well in Fig. 3, since these effects are present in both.

Finite-size scaling
The objective of finite-size scaling is to take the data in Fig. 3 and try to collapse it onto a common curve by finding suitable critical parameters. We follow the technique of ref. [25] and our previous work 1 . In particular, we use the scaling form: and we minimize a cost function with the data and its associated error as input to find the critical parameters α c, exp and ν exp . We store our experimental data as a triple α,qðαÞ,eðαÞ ð Þ , where e(α) is the standard error of the mean for the data point. The standard error of the mean for n samples is: where q i (α) is the order parameter for a given matrix B and α, and qðαÞ is the mean of q i (α) over i.
We then transform the triple according to the scaling form: We sort these triples by their t-values and then compute the cost function: with T being the number of data points. The quantity in the summation is: Δ g i À g À Á Â Ã 2 = e 2 i + t i + 1 À t i t i + 1 À t iÀ1 2 e 2 iÀ1 + t iÀ1 À t i t i + 1 À t iÀ1 2 e 2 i + 1 : The cost function Cðα c, exp , ν exp Þ measures, for each index i, the squared deviation of the point t i , g i À Á from the linear interpolation ḡ between the points t iÀ1 , g iÀ1 À Á and t i + 1 , g i + 1 À Á on either side of the sorted sequence. We exclude the first and last points in the sequence since they have no neighbouring points to the left or right, respectively. The uncertainty (Eq. (20)) is a weighted sum of the squared error of the current point t i , g i À Á and the squared error from the linear interpolation (Eq. (19)). We skip over any three identical t-values in a row in Eq. (17) because of the division by zero in Equations (19) and (20) (though this only happens for isolated values of α c, exp ). When this happens, we reduce the denominator of the fraction in front of Eq. (17) by the number of skips.
We plot the cost function over a grid of values near the critical parameters from the literature (Supplementary Fig. 2). This allows us to visualize both the minimum and the uncertainty around it. The collapse for finite-size scaling works best when there are finite-size effects, so we restricted our data for the cost function to the region α c ± 0.5 (the black connector linking the main plot with the inset in Fig. 3).
We chose our grid for the critical parameters to be α c, exp 2 0:85, 1:10 ½ , with a step size of 0.001, and ν exp 2 1:5,4:0 ½ , with a step size of 0.01. We chose a finer step size for α c, exp because we know the critical threshold. We used a larger step size for ν exp because there is less precision in the literature for ν.
To estimate our uncertainty, we plot a contour at the level 1 + r ð ÞC min , where r is the size of the maximum deviation we allow in the minimum value. We chose r = 0.25, which means we remain uncertain about the minimum for values that are up to 25% larger. Changing r will grow or shrink the contour. We note in Supplementary  Fig. 2 that the cost function's minimum resides roughly in the centre of the contour. To quantify our uncertainty, we compute the width and height of the rectangle circumscribing the contour. Then, we take the uncertainty in α c, exp to be half the width and the uncertainty in ν exp to be half the height.

Data availability
The error-mitigated output from the quantum processors is available 26 at the following Zenodo repository: https://doi.org/10.5281/zenodo. 7120441. Source data are provided with this paper.